To Eat Or Not To Eat?

Introduction

In one of the greatest cities in the world, NYC, we have an abundant amount of options to choose from when it comes to where we eat, whether we are looking for fine dining, fast food, or something in between. The cuisine choices seem just as endless. Is there a good way to help choose a restaurant to eat in or a general location? Although not an easy task, one thing we can all agree on is that nobody wants to eat in a restaurant ridden with health violations. In our project, we will explore and analyze NYC’s restaurant inspection results from 2013-2017 of the five boroughs.

We are interested in answering:

1. Which cuisines have the least and most violations along with their associated score?

2. What parts of NYC have the least and most violations?

3. We hypothesize that location is highly associated with inspection grade and so we will be searching for patterns between these variables.

4. We are also interested in seeing how inspection grade and score changes over time (years) based upon cuisine and location.

5. A description of violations is also provided in the data and we would like to better understand the common causes of those violations based upon the description.

6. Is there a pattern/trend in violations/inspection score/grade based upon restaurant chains?

We found this data set by exploring NYC open data sets. Our love for food and health made this a great option. One of our team members who recently moved to NYC is highly selective about the food she eats due to her rare health condition and hence, having the kinds of questions mentioned above answered can prove to be a great asset to both her and others suffering from health conditions or who may just be picky eaters! The data can be found and downloaded here: https://data.cityofnewyork.us/Health/DOHMH-New-York-City-Restaurant-Inspection-Results/xx67-kt59

To download, click the export button on the right -> Download as -> choose your format. For this assignement we used the CSV format.

Team

Team Members:

• Jonathan Galsurkar
• Lakshya Garg

Task Distribution

Get table format for this

library(knitr)
TaskDistribution <- read.csv(file="TaskDistribution.csv", 
                           header=TRUE, sep=",",as.is=TRUE)
kable(TaskDistribution )

S.No	Task	Assignee
1	Cleaning up the data	Jonathan Galsurkar
2	identifying missing data-including understanding what missing values mean	Lakshya Garg
3	Preparing executive summary	Jonathan Galsurkar
4	Questions 1, 2, 5	Jonathan Galsurkar
5	Questions 3, 4, 6	Lakshya Garg
6	Pushing the data analysis to GITHUB website	Both
7	Render plots in R	Both
8	Derive data for plotting in Tableau	Both
9	Creating maps in Tableau	Lakshya Garg
10	Report Writeup	Both
11	Visualization tidyup	Both

Analysis of Data Quality

To begin our analysis of data quality, let’s load our data.

ViolationsData <- read.csv(file="inspection.csv", header=TRUE, sep=",", as.is=TRUE)

Dataset Description

DatasetDescription <- read.csv(file="DescriptionDataset.csv", header=TRUE, sep=",", as.is=TRUE)
kable(DatasetDescription )

Field	Description
CAMIS	Record ID
DBA	Name of the resteraunt
BORO	Borough(Bronx,Brooklyn,Queens,Manhattan,Staten Islands)
BUILDING	Building address of resteraunt
STREET	Street address of resteraunt
ZIPCODE	ZIPCODE of resteraunt location
PHONE	Phone number fo resteraunt
CUISINE DESCRIPTION	Cusine type
INSPECTION DATE	Date of inspection-used to derive inspection year
ACTION	Action reported by inspector
VIOLATION CODE	Violation code as per blue book
VIOLATION DESCRIPTION	Violation description
CRITICAL FLAG	Criticality flag of violation reported
SCORE	Score awarded
GRADE	Grade awarder after inspection-can be missing
GRADE DATE	Date of grade award
RECORD DATE	Date the record entered the inspection results
INSPECTION TYPE	Type of inpection.initial/re-inspection

In order to more easily work with our data let’s ensure that every date is in a date format rather than a string and every string is a factor variable.

library(dplyr)
library(tidyr)
##Date conversion
ViolationsData <- ViolationsData %>%
  mutate(INSPECTION.DATE= as.Date(INSPECTION.DATE, format= "%m/%d/%Y"))%>%
  mutate(GRADE.DATE= as.Date(GRADE.DATE, format= "%m/%d/%Y"))%>%
  mutate(RECORD.DATE= as.Date(RECORD.DATE, format= "%m/%d/%Y"))
##Factor conversion
ViolationsData <- ViolationsData %>%
  mutate(CUISINE.DESCRIPTION= as.factor(CUISINE.DESCRIPTION))%>%
  mutate(BORO= as.factor(BORO))%>%
  mutate(VIOLATION.CODE= as.factor(VIOLATION.CODE))%>%
  mutate(CRITICAL.FLAG= as.factor(CRITICAL.FLAG))%>%
  mutate(GRADE= as.factor(GRADE))%>%
  mutate(INSPECTION.TYPE= as.factor(INSPECTION.TYPE))
summarize <- dplyr::summarize
mutate <- dplyr::mutate

Many of the questions we are interested in answering involve trends across restuarant locations. We first checked to see the number of restuarant inspections by borough. From the plot below, we noticed that there were a number of inspections in which the Borough information was missing and won’t help our analysis.

library(ggplot2)
boroughPlot <- ggplot(ViolationsData, aes(BORO, fill=BORO))
boroughPlot + geom_bar()+ theme(legend.position="none") +
  ggtitle('Inspection Count by Borough') + 
  labs(x = "Borough", y ="Number of Inspections")

Another main feature of our data set is the inspection year since we wish to explore patters in inspection grades/scores over the years. From the plot below, we noticed that there is almost no inspection data before 2013 and surprisingly more inspection data in 1900 than 2012 and 2011. We decided to only work with data from 2013 and up. It is important to note that the number of inspections for 2017 is low since the 2017 data is only available from January - March.

violationYearsPlot <- ggplot(ViolationsData, aes(factor(as.numeric(format(INSPECTION.DATE, '%Y'))), fill="Red"))
violationYearsPlot + geom_bar()+coord_flip()+ theme(legend.position="none") +xlab("Inspection Year")+ylab("Inspection Counts") + ggtitle("Inspection Count by Year")

Another crucial feature of our data set is grade a restaurant received after inspection. We decided to plot a stacked bar chart for to see the count of each type of grade for every cuisine. We used a stacked bar chart because we wanted to quickly assess the magnitude missing without taking up extra room. From the plot below, it was shocking to see that we generally had more missing grades than grades. This was true regardlss of cuisine. It was also interesting that the grades were purely missing and not categorized as “Not Yet Graded”.

ggplot(ViolationsData, aes(CUISINE.DESCRIPTION, fill = GRADE)) + geom_bar() + 
  coord_flip() + ggtitle("Stacked Bar Chart of Grade Distributions Across Cuisines") +
  xlab('Cuisine') + ylab('Count') +
  theme(plot.title = element_text(size = 70, face = "bold"), 
        axis.text.x = element_text(colour="grey20",size=60,angle=0,hjust=.5, 
                                   vjust=.5,face="plain"), 
        axis.text.y = element_text(colour="grey20",size=30,angle=0,hjust=1,
                                   vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=60,angle=0,hjust=.5,
                                    vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=50,angle=90,hjust=.5,
                                    vjust=.5,face="plain"),
        legend.title = element_text(size=40),
        legend.text=element_text(size=40), legend.key.size = unit(1, 'in'))

To stay on the topic of grades, after researching the letter grading program, we found the following information:

• Blank Grade (Red): The following are scored but not graded. For ex. Initial inspections that result in a score of 14 points or higher, monitoring inspections at a restaurant that has performed very poorly on its re-inspection. The Health Department may continue to inspect the restaurant roughly once a month until it scores below 28 or the Department closes it for serious and persistent violations, inspections at new restaurants not yet open to the public, an inspection at a restaurant seeking to reopen after the Department closed it, some inspections in response to complaints.

• A score of less than 14 points on either initial or re-inspection results in an “A” grade

• On re-inspection, a score of 14-27 points means a restaurant receives both a “B” grade and a “Grade Pending” card.

• On re-inspection, a score of 28 or more points means a restaurant receives both a “C” grade and a “Grade Pending” card.

• Both Z and P represent grade pending, however P represents a Grade Pending issued on re-opening following an initial inspection that resulted in a closure.

We also discovered that not every inspection is “gradable”. Gradable inspections have the following properties:

• INSPECTION TYPE in (Cycle Inspection/Initial Inspection, Cycle Inspection/Re-Inspection, Pre-Permit (Operational)/Initial Inspection, Pre-Permit (Operational)/Re-Inspection)

• ACTION in (Violations were cited in the following area(s), No violations were recorded at the time of this inspection, Establishment Closed by DOHMH)

• INSPECTION DATE > July 26, 2010

This can probably explain a fair amount of the missing grade data observed in our plot.

According to the ABOUT the data set page: The SCORE and GRADE fields may be inconsistent with each other because of limitations or errors in the data systems. That is to say, scores of 0-13, 14-27 and 28+ are not always accompanied by A, B and C grades, respectively, when they should be. There may also be cases where a grade card was given out but a record of that grade issuance is missing from the data system, and therefore missing from this dataset, even though the SCORE field is populated. Note that when initial inspections are adjudicated down to the A range, the absence of an accompanying grade associated with that inspection is correct, because the grade would not be assigned until the re-inspection is performed.

To gain some final insight on the data quality, we decided to plot the relationshp between the number of missing scores by grade and by whether or not a violation was reported. We transformed the actions taken into three categories: 1. No violations were recorded at the time of this inspection to No Violation 2. Any action reported to Violation reported 3. Missing actions to NA We then counted if the score was provided or not.

library(plyr)
score_grade <- ViolationsData[ -c(1:9, 11:13, 16:18) ]
score_grade[score_grade == ''] <- NA
score_grade_combos <- score_grade  %>% mutate(missing_score = ifelse(is.na(SCORE), "yes", "no"))
score_grade_missing <- count(score_grade_combos, c('GRADE', 'missing_score', 'ACTION'))
score_grade_missing <- score_grade_missing %>% 
  mutate(violation = ifelse(ACTION == 'No violations were recorded at the time of this inspection.', 
                                          "No violation", 
                            ifelse(is.na(score_grade_missing$ACTION), NA, "Violation Reported"))) %>%
  select(-c(ACTION))
ggplot(score_grade_missing, aes(x = GRADE, y = log(freq), fill = missing_score)) + 
  geom_bar(stat = 'identity', position = 'dodge') + facet_wrap(~violation) + 
  ggtitle("Grade and missing score combinations by violation report") +
  labs(x = "Grade", y = "Log(Frequency)")

Other interesting insights by exploring data manually are:

1. Grades of NA reported high frequency of score in both No violation and violation reported category.

2. At first, it appears that high scores are related to low grades or needs grading but then we find restaurants with a grade of A that has the same score as a restaurant with a grade of C.

3. Another insight by just looking at the data is we surprisingly saw that restaurants with a critical flag still receive grades of A.

Executive Summary

NYC puts a lot of time and money into inspecting restuarants. Why would they go through all this trouble? To ensure quality meals and satisfaction of NYC residents of course! It is important for restaurants and establishment that meet at least the minimum requirements of health and safety regulations in order to promote less food sick residents and a cleaner NYC.

How can we tell if these inspections are actually working to promote restuarants to meet regulations? Let’s take a look at the proportion of scores that were graded an A in 2013 and the proportion of scores that were graded an A in 2016. Let’s view this information by location so that you know what neighborhood to choose when you’re craving a restaurant.

The darker the shade of red, the higher proportion of A grades. It is clear that over the last 3 years, the proportion of restuarants with a grade of an A has inscreased. NYC inspections must be working to improve the quality of the restuarants we eat in! In 2013, it seems that the there were only a few neighborhoods with an extremely high proportion of A grades. Only 3 areas in Staten Island, 1 in Long Island, a few in Brooklyn, Manhattan, and the Bronx. In 2016 on the other hand, the proportion of A grades is very high almost regardless of neighborhood. Keep it up NYC restuarant inspections!

To see this from another angle, let’s look at the average scores per neighborhood in 2013 and 2016. Average scores were binned in the legend such that:

• 1 means that the average score was less than 10,

• 2 if the average score was between 10 and 20,

• 3 if the average score was between 20 and 30,

• 4 if the average score was between 30 and 40,

• 5 if it was between a 40 and 50.

Just look how much lighter the color of each neighborhood got! Lighter colors mean a lower average score, which correlates to less violations and healthier establishments!

Well now we know what neighborhoods to be in for establishments with low scores and A grades, but what are the average scores based upon cuisine?

It is interesting that Not Listed/Not Applicable and Other made it to the top 10 but there you have it, the best and worst cuisines based upon their average scores.

Main Analysis

After analyzing the quality of the data set, we got rid of data in which the Borough is missing and the year is before 2013.

ViolationsData <- ViolationsData %>% filter (BORO != "Missing")
ViolationsData <- ViolationsData %>% filter(
  as.numeric(format(INSPECTION.DATE , '%Y')) > 2012)

Now we can get a better picture of the total number of inspections by year and borough.

boroughPlot <- ggplot(ViolationsData, aes(BORO,fill=BORO))
boroughPlot + geom_bar() + theme(legend.position="none") + 
  facet_wrap(~factor(as.numeric( format(INSPECTION.DATE , '%Y')))) + 
  ggtitle("Total Inspections from 2013-2017 in each Borough") +
  xlab("Borough") + ylab("Number of Inspections") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

After 2013, there seems to be a consistent amount of inspections across the years. The number of inspections by borough also seems to make sense since we expect Manhattan to have the largest number of restuarants.

Next, we took a look at the grade distribution by borough.

inspection_grades <- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, -ZIPCODE, -PHONE, -ACTION, -VIOLATION.CODE, -VIOLATION.DESCRIPTION, -CRITICAL.FLAG, -SCORE, -GRADE.DATE, -RECORD.DATE, -INSPECTION.TYPE,-CUISINE.DESCRIPTION)
inspection_grades_woyear<-inspection_grades %>% select(-INSPECTION.DATE)
inspection_grades_woyear <- inspection_grades_woyear %>% gather(key, value, -BORO) %>% group_by(BORO, key, value) %>% tally %>% spread(value, n, fill = 0)
names(inspection_grades_woyear)[names(inspection_grades_woyear)=="key"] <- "grade"
names(inspection_grades_woyear)[names(inspection_grades_woyear)==""] <- "`Unknown"
inspection_grades_woyear <- inspection_grades_woyear %>% gather(key, value, -BORO,-grade)
RelFreq<-function(m){
   ((m )/sum(m))
}
inspection_grades_woyear<-inspection_grades_woyear %>% group_by(BORO,grade) %>%
    mutate(percentage = (value/sum(value))*100)
ggplot(inspection_grades_woyear, aes(BORO, percentage ,fill= key)) + 
  geom_bar(stat="identity", position = "dodge") +
  ggtitle('%Grade Distribution across boroughs') + 
  labs(x = "Borough", y ="Percentage") + guides(fill=guide_legend(title="Grade"))  + theme(plot.title = element_text(size = 25, face = "bold"))

We initially hypothesized that we would see a pattern in grade distribution by borough, however, the plot shows us that grade distributions are almost the same regardless of borough.

inspection_grades_year<-inspection_grades%>% mutate(year=factor(as.numeric(format(INSPECTION.DATE , '%Y'))))
inspection_grades_year<-inspection_grades_year%>% select(-INSPECTION.DATE)
 inspection_grades_year <- inspection_grades_year %>% gather(key, value, -BORO,-year) %>% group_by(BORO,year, key, value) %>% 
   tally %>% spread(value, n, fill = 0)# %>% gather(blah, -BORO, -key) #summarize(aprop = A/(A+B+C+))
names(inspection_grades_year)[names(inspection_grades_year)=="key"] <- "grade"
names(inspection_grades_year)[names(inspection_grades_year)==""] <- "`Unknown"
inspection_grades_year <- inspection_grades_year %>% gather(key, value, -BORO,-grade,-year)
RelFreq<-function(m){
   ((m )/sum(m))
 }
inspection_grades_year<-inspection_grades_year %>% group_by(BORO,grade,year) %>%
    mutate(percentage = (value/sum(value))*100)
ggplot(inspection_grades_year, aes(BORO, percentage ,fill= key)) + 
  geom_bar(stat="identity", position = "dodge") +facet_wrap(~year,nrow=2,ncol=3)+
  guides(fill=guide_legend(title="Grade"))+
  ggtitle("%Grade Distribution by Borough and Year")+ theme(plot.title = element_text(size = 25, face = "bold"))

Adding year to this analysis showed us an increase in the proportian of As in Staten Island in 2015 but seemed consistent throughout the rest of the plot.

We decided that borough may be too general and thus looked at grade distribution by zip code across the various years. We used a heat map to do so.

library(viridis)
nonYearDataForHeatMap<-ViolationsData %>% select(GRADE,ZIPCODE)
nonYearDataForHeatMap <- nonYearDataForHeatMap %>% group_by(GRADE,ZIPCODE) %>% tally 
RelFreq<-function(m){
   ((m )/sum(m))
}
nonYearDataForHeatMap<-  nonYearDataForHeatMap %>% group_by(ZIPCODE)%>%
    mutate(percentage = (n/sum(n))*100)
ggplot(nonYearDataForHeatMap, aes(GRADE, 
                        ZIPCODE, fill = percentage)) +
  geom_tile() +
  scale_fill_viridis() +
  ggtitle("Percentage Grade Distribution by Zipcode\n ")+ylab("Zip Code") + theme(plot.title = element_text(size = 35), 
        axis.text.x = element_text(colour="grey20",size=15,angle=0,hjust=.5, 
                                   vjust=.5,face="plain"), 
        axis.text.y = element_text(colour="grey20",size=10,angle=0,hjust=1,
                                   vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=25,angle=0,hjust=.5,
                                    vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=30,angle=90,hjust=.5,
                                    vjust=.5,face="plain"),
        legend.title = element_text(size=20),
        legend.text=element_text(size=20), legend.key.size = unit(1, 'in'))

It is interesting to note that missing grades are probably the most probable in general. There are some zip codes with almost no missing grades. The restuarants in those same zip codes seem to mostly have all A grades. Another interesting point is 1 zip code which has all Not Yet Graded grades. In general missing grades and grades of an A seem the most common. This is followed by grades of a B.

The next part of our analysis was to look at the average scores of each cuisine. With this information, we can help consumers see what kinds of establishments have the best and worst scores on average. This in term can help a consumer choose a type of cuisine when they are hungry.

average_scores <- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, -ZIPCODE, 
                                            -PHONE, -INSPECTION.DATE, -ACTION, 
                                            -VIOLATION.CODE, -VIOLATION.DESCRIPTION, 
                                            -CRITICAL.FLAG, -GRADE, -GRADE.DATE, 
                                            -RECORD.DATE, -INSPECTION.TYPE)
average_cuisine <- average_scores %>%  group_by(CUISINE.DESCRIPTION) %>% na.omit %>%
  dplyr::summarize(average_score = mean(SCORE))
average_cuisine <- arrange(average_cuisine, -average_score)
ggplot(average_cuisine, aes(reorder(x = CUISINE.DESCRIPTION, --average_score), average_score,fill="Blue")) + geom_bar(stat='identity') + coord_flip() + theme(legend.position = 'none')+ggtitle("Average Scores Across Cuisines")+xlab("Cuisine Description")+ylab("Average Score")+geom_text(aes(label = sprintf("%.2f", average_score)),position=position_stack(vjust=0.5),vjust = 0.5, size = 3) + theme(plot.title = element_text(size = 25), 
        axis.title.x = element_text(colour="grey20",size=20,angle=0,hjust=.5,
                                    vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=20,angle=90,hjust=.5,
                                    vjust=.5,face="plain"))

worst10 <- average_cuisine[1:10,]
ggplot(worst10, aes(reorder(x = CUISINE.DESCRIPTION, --average_score), average_score, fill = CUISINE.DESCRIPTION)) + geom_bar(stat='identity') + coord_flip() + theme(legend.position = 'none')+ggtitle("Average Score of Worst 10 cuisines")+xlab("Cuisine Description")+ylab("Average Score")+geom_text(aes(label = sprintf("%.2f", average_score)),position=position_stack(vjust=0.5))

n <- nrow(average_cuisine)
best10 <- average_cuisine[(n- 10):n,]
best10 <- best10[order(best10$average_score),]
ggplot(best10, aes(reorder(x = CUISINE.DESCRIPTION, -average_score), average_score, fill = CUISINE.DESCRIPTION)) + geom_bar(stat='identity') + coord_flip() + theme(legend.position = 'none')+ggtitle("Average Score of top 10 cuisines")+xlab("Cuisine Description")+ylab("Average Score")+geom_text(aes(label = sprintf("%.2f", average_score)),position = position_dodge(width = 1),#position=position_stack(vjust=0.5)
                                                                                                                     vjust = 0.5, size = 4)

It is interesting that the mysterious “not applicable”" cuisine is one of the best!

In this next plot, we wanted to see if the average score significantly differed by Borough.

average_borough <- average_scores %>%  group_by(BORO) %>% na.omit %>%
  dplyr::summarize(average_score = mean(SCORE))
ggplot(average_borough, aes(reorder(x = BORO, --average_score), average_score, fill = BORO)) + geom_bar(stat='identity') + theme(legend.position = 'none')+ggtitle("Average Score across boroughs")+xlab("Borough")+ylab("Average Score")+geom_text(aes(label = sprintf("%.2f", average_score)),position = position_dodge(width = 1),
                                                                                                                     vjust = -0.5, size = 3)

Although we see a slight difference, the boroughs have approximately the same average score. Staten Island may be the one exception. We also faceted by year to see if that made a difference, but the results were essentially the same as in the plot above.

So far we looked at grades and scores but not the actual violations. Whare are the top violations restaurants usually face?

violations<- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, -ZIPCODE, 
                                            -PHONE, -INSPECTION.DATE, -ACTION, 
                                            -VIOLATION.CODE, 
                                            -CRITICAL.FLAG, -GRADE, -GRADE.DATE, 
                                            -RECORD.DATE, -INSPECTION.TYPE)
vs <- violations %>%  group_by(VIOLATION.DESCRIPTION) %>% na.omit %>%
  dplyr::summarize(count = n())
vs <- arrange(vs, -count)
topviolations <- vs[1:10,]
library(stringr)
topviolations$viol = str_wrap(topviolations$VIOLATION.DESCRIPTION, width = 15)
ggplot(topviolations, aes(reorder(viol, -count), count,fill=viol)) + geom_bar(stat='identity') + 
  theme(legend.position = 'none')+ggtitle("Top 10 violations \n")+xlab("Violation description")+ylab("count")+geom_text(aes(label = sprintf("%.0f", count)),
 position = position_dodge(width = 1), vjust = -0.5, size = 7) + theme(plot.title = element_text(size = 40, face = "bold"), 
        axis.text.x = element_text(colour="grey20",size=20,angle=0,hjust=.5, 
                                   vjust=.5,face="plain"), 
        axis.text.y = element_text(colour="grey20",size=30,angle=0,hjust=1,
                                   vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=30,angle=0,hjust=.5,
                                    vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=30,angle=90,hjust=.5,
                                    vjust=.5,face="plain"))

We also hypothesized that the type of inspections would generally vary by Borough but even this was mostly consistent.

 library(viridis)
mos<- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, -ZIPCODE, 
                                           -PHONE, -ACTION, 
                                           -VIOLATION.CODE, 
                                           -CRITICAL.FLAG, -GRADE, -GRADE.DATE, 
                                           -RECORD.DATE, -SCORE, -CUISINE.DESCRIPTION, -VIOLATION.DESCRIPTION)
mos$year <- factor(as.numeric( format(mos$INSPECTION.DATE , '%Y')))
mos <- mos %>% select(-INSPECTION.DATE)
average_mos <- mos %>% group_by(BORO, year, INSPECTION.TYPE) %>% tally %>%
 group_by(BORO, year)  %>% 
 mutate(percentage = n / sum(n))
#average_cuisine <- arrange(average_cuisine, -average_score)
library(vcd)
average_mos <- average_mos %>% select(-n)
ggplot(average_mos, aes(BORO, 
                       INSPECTION.TYPE, fill = percentage)) +
 geom_tile() +
 scale_fill_viridis() +
 #facet_wrap(~BORO)
 ggtitle("Average Violation Score by Zip Code\n ")+ xlab('Borough') +
 theme(plot.title = element_text(size = 30, face = "bold"), 
        axis.text.x = element_text(colour="grey20",size=20,angle=0,hjust=.5, 
                                   vjust=.5,face="plain"), 
        axis.text.y = element_text(colour="grey20",size=15,angle=0,hjust=1,
                                   vjust=0,face="plain"),  
        axis.title.x = element_text(colour="grey20",size=25,angle=0,hjust=.5,
                                    vjust=0,face="plain"),
        axis.title.y = element_text(colour="grey20",size=25,angle=90,hjust=.5,
                                    vjust=.5,face="plain"),
        legend.title = element_text(size=30),
        legend.text=element_text(size=20), legend.key.size = unit(1, 'in'))

One thing we noticed was that Staten island has the least Cycle Inspection/Initial Inspection. Maybe they were opening up the least new restuarants. They also had the most cycle inspection/re-inspection percentage. As we’ll see in the maps below, those reinspections may have helped restaurants in staten island improve their grades and scores.

One factor of our data set we haven’t discussed too much is the Critical flag, saying whether or not a violation is critical. Once again, we assumed that over the years and across different boroughs, there would be some pattern in critical violations. As can be inferred from the following plot, there is not.

CriticalityData<- ViolationsData %>% select(BORO,INSPECTION.DATE,CRITICAL.FLAG) %>% na.omit()
CriticalityData$year<-factor(as.numeric( format(CriticalityData$INSPECTION.DATE , '%Y')))
#CriticalityData <- CriticalityData %>%select(-INSPECTION.DATE)
CriticalityData <- CriticalityData %>%select(-INSPECTION.DATE) %>% gather(key, value, -BORO, -year) %>% group_by(BORO, year,key,value) %>% tally %>% spread(value, n, fill = 0)
names(CriticalityData)[5] <- "NotApplicable"
names(CriticalityData)[6] <- "NonCritical"
CriticalityDataPercentage<-CriticalityData %>% summarize(criticalPercent = Critical/(Critical+NotApplicable+NonCritical), NotApplicablePercent = NotApplicable/(Critical+NotApplicable+NonCritical), NonCriticalPercent = NonCritical/(Critical+NotApplicable+NonCritical))#
names(CriticalityDataPercentage)[4] <- "Critical"
names(CriticalityDataPercentage)[5] <- "NotApplicable"
names(CriticalityDataPercentage)[6] <- "NonCritical"
CriticalityDataPercentage<-CriticalityDataPercentage%>% gather(key1, value, -BORO, -year,-key)
ggplot(CriticalityDataPercentage, aes(BORO, value, fill = key1)) + 
  geom_bar(stat="identity", position = "dodge") + 
  scale_fill_brewer(palette = "Set1")+ xlab("BOROUGH")+ylab("CRITICALITY PERCENTAGES") +theme(axis.text.x = element_text(angle = 45, hjust = 1))+ ggtitle("Grouped barchart of grade frequency in each cuisine category \n when there was a violation reported")+facet_wrap(~year,nrow=2,ncol=3)+guides(fill=guide_legend(title="Criticality \n category"))

We attempted to further analyze the critical flag but we learned that the critical flag did not have too much of an impact on grade or scores. There was a lot of data with a grade of A but a Critical flag while others with a very high score by no critical flag.

VISUALISING DATA WITH MAPS

For the remainder of the analysis, we will be working with Tableau. We chose Tableau because we wanted to plot elegent maps. Tableau provides a simpler and more elegent way to do so than in R.

• Please note that there is a tutorial for deriving maps from generated datset in Tableau in the github repository as ChoroplethMaps.html under the “FinalReport” folder in the Github Repository. It was created by our very own Lakshya Garg.

The code chunks will generate the data that we upload to Tableau. In the following analysis, we explore the proportion of A grades based on neighborhoods over the years.

#####GENERATE GRADE PROPOERTION DATASET BY YEAR FOR TABLEAU PLOTTING
zipdata<- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, 
                                            -PHONE, -ACTION, 
                                            -VIOLATION.CODE, 
                                            -CRITICAL.FLAG, -SCORE, -GRADE.DATE, 
                                            -RECORD.DATE, -INSPECTION.TYPE, -VIOLATION.DESCRIPTION, -CUISINE.DESCRIPTION, -BORO) %>% na.omit()
zipdata$year<-factor(as.numeric( format(zipdata$INSPECTION.DATE , '%Y')))
zipdataGradePercentage <- zipdata %>%select(-INSPECTION.DATE) %>%gather(key, value, -ZIPCODE, -year) %>% group_by(ZIPCODE, year, key, value) %>%
  tally %>% spread(value, n, fill = 0) %>% summarize(aprop = A/(A+B+C), bprop = B/(A+B+C), cprop = C/(A+B+C) )
# ggplot(zipdataGradePercentage, aes(reorder(x = ZIPCODE, --aprop), aprop)) + geom_bar(stat='identity') + 
#   coord_flip() + facet_wrap(~year)
#   theme( axis.text = element_text(size = 3))

PROPORTIONS OF A GRADES ACROSS THE YEARS

Over the years, the proportion of A’s has increased, as represented by the darker red shading of the map. It is too early to study 2017 in depth due to the limited number of inspections so far. Note that 2017 data maps are only for the first quarter and are likely to change.

knitr::include_graphics('Aprop2013.png')

knitr::include_graphics('Aprop2014.png')

knitr::include_graphics('Aprop2015.png')

knitr::include_graphics('Aprop2016.png')

knitr::include_graphics('Aprop2017.png')

PROPORTIONS OF B GRADES ACROSS THE YEARS

Let’s look at how the proportion of B’s change. We can see that over the years, there seem to be less B’s which makes sense since over the years we have more grades of an A. Note that 2017 data maps are only for the first quarter and are likely to change.

knitr::include_graphics('Bprop1.png')

knitr::include_graphics('Bprop2.png')

knitr::include_graphics('Bprop3.png')

knitr::include_graphics('Bprop4.png')

knitr::include_graphics('Bprop5.png')

AVERAGE SCORES BY ZIPCODE OVER THE YEARS

Let’s look at the average scores by neighborhood. We have already done so by cuisine. Over the years, we have a lower Average Score. These NYC inspections must be prompting resturants to improve their facilities and follow regulations!

Note that 2017 data maps are only for the first quarter and are likely to change.

knitr::include_graphics('AvgScore2013.png')

knitr::include_graphics('AvgScore2014.png')

knitr::include_graphics('AvgScore2015.png')

knitr::include_graphics('AvgScore2016.png')

knitr::include_graphics('AvgScore2017.png')

AVERAGE RE-INSPECTIONS OVER THE YEARS

We decided to see which areas had the highest average number of reinspections. Maybe this can explain the better scores and grades over time. Interestingly enough, the areas with the most reinspections seem to be the ones in which the grades and scores improved. The system of reinspections must be working!

Note that 2017 data maps are only for the first quarter and are likely to change.

### DATA GENERATION
library(stringr)
reinspect <- ViolationsData %>% select(-CAMIS, -DBA, -BUILDING, -STREET, -CUISINE.DESCRIPTION, -BORO, -SCORE,
                                           -PHONE, INSPECTION.DATE, -ACTION, 
                                           -VIOLATION.CODE, 
                                           -CRITICAL.FLAG, -GRADE, -GRADE.DATE, 
                                           -RECORD.DATE)
reinspect$year<-factor(as.numeric( format(reinspect$INSPECTION.DATE , '%Y')))
#inspection_types <- INSPECTION.TYPE
reinspect <- reinspect %>%select(-INSPECTION.DATE)  %>% na.omit %>% filter(str_detect(INSPECTION.TYPE, "Re-inspection"))%>%group_by(VIOLATION.DESCRIPTION, year, ZIPCODE )  %>% tally# summarise(re = sum(str_count(reinspect$INSPECTION.TYPE, "Re-inspection")))

knitr::include_graphics('AvgReinspect2013.png')

knitr::include_graphics('AvgReinspect2014.png')

knitr::include_graphics('AvgReinspect2015.png')

knitr::include_graphics('AvgReinspect2016.png')

knitr::include_graphics('AvgReinspect2017.png')

STARBUCKS OR DUNKIN DONUTS???

As a last bit of analysis, we decided to help settle a huge dilemna : Where should you get your coffee, Starbucks or Dunkin Donuts?

##DATA GENERATION CODE
library(stringr)
ChainssDf<-ViolationsData%>% mutate(isStarbucks=str_detect(DBA,"STARBUCK")) %>% mutate(isDunkin=str_detect(DBA,"DUNKIN"))
ChainssDf$DBA[ChainssDf$isStarbucks==TRUE]<-"STARBUCKS"
ChainssDf$DBA[ChainssDf$isDunkin==TRUE]<-"DUNKIN"
ChainssDf<-ChainssDf%>%filter(DBA %in% c("STARBUCKS","DUNKIN"))
ChainssDf<-ChainssDf%>%select(DBA,BORO,INSPECTION.DATE,ZIPCODE,SCORE,VIOLATION.DESCRIPTION,GRADE)
ChainssDf$year<-factor(as.numeric( format(ChainssDf$INSPECTION.DATE , '%Y')))
ChainssDf<-ChainssDf%>%select(-INSPECTION.DATE)

As you can see, Dunkin Donuts has more critical violations than Starbucks in every Borough. Noticible differences are especially see in Brooklyn and Manhattan.

knitr::include_graphics('DS1.png')

Both establishments shared the same top violations. Let’s see what percentage of these violations each has.

knitr::include_graphics('DS2.png')

Once again, Starbucks is the winner!

What about the average score?

knitr::include_graphics('DS3.png')

Regardles of the year, the average score of Dunkin Donuts has been higher. The higher the score, the worse. Noticiable differences are seen in 2015 and the first three months of 2017.

Let’s see if the Borough you’re in should impact your choice.

knitr::include_graphics('DS4.png')

We guess when you’re in Brooklyn, it doesn’t matter but everywhere else, stick to Starbucks, especially you Staten Island folks!

Conclusion

Many limitations included the missing data itself. We learned that sometimes, what may seem like an intuitive hypothesis, is actually completely false. There were not as many patterns as we anticipated in terms of Borough. The only thing we did see is that over time, the score and grade has increased, showing us that the inspections seem to be working and that restuarants are improving their facilities to receieve a better score and grade. According to the website we received the data from, three fields are soon to be added as data logic becomes available to populate them accurately. Those fields are VIOLATION POINTS (the points assigned to a violation before or after adjudication, depending on whether adjudication has occurred), FINES TOTAL (the fine amount after adjudication), and DECISION DATE (adjudication date – or date a grade becomes final). With this information, we can do more analysis to see the fines of various violations and learn how each violation actually affects the score. With that, restaurants can focus on the main violations to improve their scores in the next inspection.